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1.  INTRODUCTION 


Accurate  knowledge  of  electron  densities  in  the  ionosphere  is  important  for  enabling  fast  high-accuracy 
positioning,  reliable  ground-to-satellite  communications,  and  surveillance.  Thus,  an  ability  to  forecast 
global  and  regional  ionospheric  conditions  has  applications  in  many  areas  of  military  operations,  from 
HF  communications  for  delivering  vital  information  to  soldiers  in  remote  battlefields,  to  precision- 
guided  weapons,  space-based  intelligence  gathering,  and  space-object  tracking. 

The  atmosphere,  including  the  ionosphere,  is  a  chaotic  system;  small  errors  in  the  initial  conditions  of  a 
forecast  can,  under  some  circumstances,  grow  rapidly  and  affect  predictability.  Furthermore, 
predictability  is  limited  by  errors  due  to  the  approximate  simulation  of  relevant  physical  processes  in 
the  numerical  models  and  to  poorly  known  external  forcing.  A  process  known  as  data  assimilation  aims 
to  decrease  these  uncertainties  by  using  observations  in  order  to  obtain  better  initial  conditions  and/or  to 
provide  better  estimates  of  poorly  known  empirical  quantities  in  parameterizations  of  various  physical 
processes  in  the  models.  Fundamentals  of  data  assimilation  were  developed  for  numerical  weather 
prediction  (NWP)  and  have  been  in  use  for  some  time.  Modem  NWP  models  include  1-10  million  grid- 
points  and  initial  atmospheric  states  (temperature,  humidity,  etc)  need  to  be  specified  at  each  of  these 
grid-points  for  forecasting.  It  is  clearly  impossible  to  have  an  observational  network  of  this  size  co¬ 
located  with  model  grid-points.  Additionally,  meteorological  data  (both  ground-based  and  satellite- 
collected)  often  have  large  and  sometimes  poorly  known  systematic  and  random  errors. 

Significant  improvements  in  modem  NWP  have  been  driven,  in  large  part  by  advances  in  the  data 
assimilation  methodology.  In  very  simplistic  terms,  data  assimilation  aims  to:  (1)  use  model- 
encapsulated  knowledge  of  underlying  atmospheric  physics  to  interpolate  between  observations;  (2) 
use  Bayesian  statistics  to  continuously  combine  model  state  with  observations  thus  producing  a  less 
uncertain  state;  and  (3)  compute  temporal  evolution  and  spatial  distribution  of  model  errors  (error 
covariances)  to  keep  track  of  forecast  quality  and  to  be  able  to  reject  (or  correct)  inconsistent 
observations. 

In  order  to  apply  this  methodology  to  ionospheric  forecasting  one  needs  a  numerical  model  of  the 
ionosphere.  While  present  meteorological  forecast  models  are  very  mature,  driven  by  the  need  for  a 
reliable  weather  forecasting,  ionospheric  numerical  models  are  mainly  academic  and  not  as  advanced. 
Even  less  advanced  are  the  specific  data  assimilation  techniques  that  can  be  applied  to  the  ionosphere. 
Additionally,  in  meteorology  one  has  to  solve  fluid  dynamics  equations  for  only  two  constituents,  air 
and  water  vapor.  Modem  ionospheric  physics  has  to  consider  separately  seven  different  ions,  all  free 
electrons,  seven  neutral  particles,  and  the  interactions  between  them  which  adds  significant  complexity. 

The  amount  of  ions  and  electrons  in  a  particular  region  of  the  ionosphere  is  governed  by  a  multitude  of 
complex  physical  and  chemical  processes.  Charged  particles  appear  when  Sun-emitted  photons  strike 
neutral  atoms  or  molecules  in  the  atmosphere  and  ionize  them.  Since  the  sunlight  intensity  increases 
with  altitude  but  the  amount  of  neutral  particles  decreases,  electron  concentrations  peak  at  about  300 
km.  Once  ions  and  electrons  have  appeared,  they  interact  with  each  other  and  neutral  particles  along 
with  the  Earth’s  magnetic  field  in  a  complex  way.  Additional  processes  need  to  be  considered  in  the 
polar  and  equatorial  regions.  All  these  interactions  result  in  fast  temporal  changes  and  often  sharp 
horizontal  gradients  in  the  electron  content;  an  excellent  review  of  ionospheric  processes  can  be  found 
in  Schunk  and  Nagy,  (2000).  Building  a  numerical  model  simulating  all  these  processes  and  assuring 
that  the  resulting  code  is  computationally  efficient  is  a  difficult  task. 

This  report  describes  a  global  three-dimensional  numerical  model  of  the  ionosphere  and  a  data 
assimilation  methodology  developed  by  Fusion  Numerics  Inc.  The  model,  sponsored  by  the  US  Air 
Force,  addresses  GPS  receiver  data  assimilation  and  bias  estimation  methodologies.  A  part  of  the 
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system  has  been  operational  since  August  2003  and  can  provide  real-time  slant  TEC  from  any  location 
on  Earth  to  all  visible  GPS  satellites:  http://fusionnumerics.com 'ionosphere. 

Similar  research  projects  have  been  undertaken  earlier  as  a  part  of  the  Air  Force  Multi-Disciplinary 
University  Research  Initiative  GAIM  (Global  Assimilation  of  Ionospheric  Measurements;  e.g., 

http://genesis2.ipl.nasa.gov/archive/2002 1 203 1  , 

http://m\As\cosismet/abstracts/EAE03/ 1 26  i 4EAE03-J- 1 26 1 4.pdf).  The  current  project  is  somewhat 
parallel  to  that  effort.  We  addressed  quite  differently  a  number  of  important  issues  related  to  numerical 
representation  of  ionospheric  physics  in  the  model  and  implementation  of  data  assimilation.  Some  of 
our  main  objectives  were  creating  a  modem,  well  documented,  and  maintainable  numerical  ionospheric 
model  and  its  early  operational  implementation.  Note  that  as  of  this  writing  ionospheric  specifications 
for  the  Air  Force  are  provided  operationally  by  PRISM  (Parameterized  Real-time  IonoSpheric  Model). 
One  of  die  objectives  of  the  Air  Force  program  solicitations  that  funded  this  effort  was  to  build  a 
prototype  of  the  next  generation  operational  ionospheric  model  capable  of  assimilating  different  types 
of  real-time  data. 

Implementation  of  practical  and  efficient  data  assimilation  schemes  as  well  as  operational 
implementations  of  forecast  systems  depends  critically  on  the  design,  quality,  and  maintainability  of  the 
underlying  physical  propagator  model.  During  the  first  year  of  this  AFRL-sponsored  investigation  we 
developed  a  new  computer  model  for  simulating  time  evolution  of  ion  and  electron  densities  in  the 
ionosphere  on  a  global  scale  and  developed  the  corresponding  tangent  linear  and  adjoint  models  with 
respect  to  parallel  transport.  The  core  ionospheric  model  solves  plasma  dynamics  and  composition 
equations  governing  evolution  of  density,  velocity  and  temperature  for  seven  ion  species  and  electrons 
on  a  fixed  grid  in  magnetic  coordinates.  It  uses  a  realistic  model  of  the  Earth's  magnetic  field  and  solar 
indices  obtained  in  real  time  from  NOAA’s  Space  Environment  Center.  At  the  present  time  the  model 
computes  real-time  ion  and  electron  densities  at  a  grid  of  more  than  one  million  points.  Higher 
resolutions  are  anticipated  in  the  future.  During  the  second  year,  we  developed  and  tested  a  new  large- 
scale  suboptimal  Kalman  Filter  data  assimilation  module  for  systematically  adjusting  a  model-derived 
state  (3-D  fields  of  electron  densities)  with  slant  total  electron  content  (TEC)  measured  by  a  network  of 
dual  frequency  geographically  distributed  GPS  receivers. 

Results  of  the  2-year  investigation  are  described  in  this  report.  Some  of  these  results  have  also  been 
presented  at  the  IEEE  Position,  Location  and  Navigation  2004  (PLANS  2004)  symposium  (Khattatov 
et  al,  2004),  and  European  Navigation  Conference  2004  (ENC  2004).  A  published  paper  from  the 
proceedings  of  PLANS  2004  can  be  found  in  the  appendix  of  this  report. 

In  the  beginning  of  this  investigation  we  decided  that  the  computer  code  of  the  ionospheric  model 
available  to  us.  Coupled  Thermosphere  Ionosphere  Model  (CTIM),  was  unacceptable  for  the  purpose 
of  adding  data  assimilation  capabilities  and  practical  implementation  of  the  forecast  system  for  several 
reasons,  with  the  major  reasons  listed  below: 

Very  poor  quality  of  code  leading  to  inefficiencies  and  redundancy,  rendering  the  code  hard  to 
understand,  maintain,  and  modify. 

Almost  complete  lack  of  documentation. 

Extensive  use  of  FORTRAN  77  and  FORTRAN  66  features  that  are  being  phased  out  in  the 
modem  and  future  FORTRAN  compiler  versions. 

Thus,  the  first  major  task  of  this  investigation  was  creating  a  new  set  of  computer  code  implementing 
the  forward  model  and  its  validation.  Certain  features  of  the  model  that  are  deemed  as  non-critical 
remain  in  the  development  stage  (e.g.,  correct  calculations  of  night  time  ionization  rates)  and  will  be 
added  later.  While  development  of  the  propagator  model  “from  scratch”  took  time,  we  deem  it  to  be  a 
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critical  component  that  enabled  and  strongly  facilitated  reaching  the  final  objective  -  building  a  modem 
and  operational  data  assimilation  and  forecasting  system  of  the  ionosphere. 

Traditionally,  development  of  such  models  in  academia  took  decades  and  the  efforts  of  many  scientists. 
While  resulting  computer  codes  would  be  cutting  edge  in  the  scientific  sense,  they  usually  would  be  not 
particularly  computationally  efficient,  user  friendly  or  well  documented.  Consequentially,  it  would  be 
difficult  to  maintain  and  advance  such  codes  or  build  operational  systems  based  upon  them.  Our 
approach  involved  reviewing  key  scientific  publications  describing  relevant  ionospheric  physics  and 
capturing  the  described  concepts  and  equations  in  a  series  of  “use-cases.”  These  use-cases  were  later 
coded  by  professional  software  engineers  and  computer  scientists  in  a  modem  object-oriented  fashion. 

As  a  result  of  this  approach,  which  clearly  would  not  be  possible  without  many  years  of  government- 
sponsored  academic  research,  energy,  and  efforts  of  many  devoted  scientists  around  the  world,  Fusion 
Numerics  has  built  a  modem  numerical  ionospheric  model  in  a  little  over  a  year  and  then  coupled  it 
with  an  efficient  data  assimilation  scheme.  Web-based  access  to  the  system  is  provided  to  early  users 
for  validation  and  exploration  purposes  at  http://fusionnumerics.com/ionosphere. 


2.  METHODS,  ASSUMPTIONS,  AND  PROCEDURES 

2.1  Model 

The  developed  model  is  a  numerical  global  model  of  the  ionosphere  loosely  based  on  a  description 
given  in  Millward  et  al  (1996)  and  Bailey  and  Balan  et  al  (1996).  It  consists  of  over  60,000  lines  of 
C++  code  and  was  ported  to  SuSe  Enterprise  64  bit  Linux,  SuSe  and  RedHat  32  bit  Linux  as  well  as 
MS  Windows  platforms.  Since  the  code  is  fully  64-bit  compliant  and  adheres  to  standard  C++  language 
specifications  it  is  expected  to  run  on  all  major  UNIX  systems. 

The  code  has  been  developed  using  modem  software  engineering  tools  and  methodologies  (OO 
design,  Rational  Unified  Process,  unit  testing,  UML  diagramming  and  reverse  engineering)  and 
profiled  and  optimized  with  industry  standard  profilers  (Rational  Rose  and  Intel  VTune).  While 
adherence  to  industry  standard  large-scale  software  development  process  allowed  us  to  develop  the 
model  in  a  relatively  short  time  frame,  this  project  would  not  have  been  possible  without  decades  of 
published  academic  research  by  many  scientists  worldwide  and  respective  financial  support  from 
various  sponsors. 

The  model  computes  spatial  distribution  and  temporal  evolution  of  seven  major  ions  (H+,  0+,  He+, 
O2+,  NO+,  N2+,  N+)  and  electrons.  Other  prognostic  variables  include  ion  and  electron  temperatures 
and  velocities.  The  model  solves  plasma  dynamics  equations  for  all  seven  ion  species  and  electrons  and 
an  energy  conservation  equation  for  the  three  major  ions  and  electrons.  It  includes  chemical  interactions 
with  neutrals  and  ions,  recombination,  ion-ion  and  ion-neutral  collision  rates,  photoionization,  and 
different  types  of  heating.  The  model  variables  are  described  in  Table  1,  Table  2  and  Table  3. 

The  model  domain  covers  all  latitudes  and  longitudes;  however  implementation  of  polar  transport  and 
high-latitude  effects  terms  are  still  in  the  experimental  stage  and  should  not  be  relied  upon.  As 
customary  in  ionospheric  applications,  the  dynamic  equations  are  solved  in  so-called  magnetic 
coordinates  since  in  the  absence  of  electric  fields  plasma  predominantly  moves  parallel  to  the  direction 
of  the  magnetic  field.  A  detailed  discussion  of  the  coordinate  transformation  and  related  equations  can 
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be  found,  for  instance,  in  Millward  et  al  ( 1 996)  and  Bailey  and  Balan  et  al  ( 1 996).  Here  we  give  only  a 
brief  overview  for  the  benefit  of  readers  not  familiar  with  this  subject. 

The  Earth’s  magnetic  field  is  being  approximated  as  that  of  a  tilted  eccentric  dipole.  The  first  magnetic 
coordinate  is  magnetic  longitude.  For  each  magnetic  longitude  we  consider  a  “stack”  of  magnetic  field 
lines  characterized  by  the  distance  of  the  apex  of  each  line  from  the  Earth’s  center  at  the  magnetic 
equator.  This  distance,  normalized  by  the  Earth’s  radius,  is  the  second  magnetic  coordinate,/?.  Finally, 
for  each  field  line  the  distance  from  the  apex  to  a  particular  point  along  the  line  gives  the  third 
coordinate,  s. 


A  portion  of  the  model  grid  for  low-latitudes  only  is  shown  in 
Figure  1.  For  clarity  only  20  longitudes  and  30  p  values  are  shown. 

The  model  solves  several  types  of  prognostic  equations.  These  equations  are  given  in  dipole 
coordinates,  along  magnetic  field  lines.  Therefore  there  is  only  one  dependent  spatial  coordinate 
corresponding  to  the  position  along  the  magnetic  field  line,  s. 

Once  the  field-aligned  transport  is  solved,  we  compute  plasma  evolution  due  to  cross-field  transport. 
Cross-field  transport  is  forced  by  electric  fields  either  imposed  externally  from  the  magnetosphere  or 
generated  internally  from  the  action  of  the  neutral  wind.  In  the  lower  thermosphere  the  mobility  of  the 
ions  is  inhibited  by  collisions  with  the  neutral  atmosphere.  The  dynamo  action  of  the  neutral  winds 
drives  currents  that  through  continuity  create  polarization  electric  fields.  The  ions  respond  to  these 
electric  fields  by  drifting  perpendicular  to  both  the  electric  and  magnetic  fields. 

This  is  often  referred  to  as  ExB  transport.  In  non-fully  coupled  models  the  ExB  plasma  velocity  is 
specified  from  external  empirical  models.  Once  this  velocity  is  known,  solving  plasma  advection 
equation  is  relatively  straightforward.  To  model  ExB  drif  at  low  latitudes  we  use  coefficients  of  the 
Fejer  and  Schierless  model.  To  model  ExB  drift  in  the  polar  regions  we  use  coefficients  of  the  Weimer 
model  provided  by  Dr.  Dan  Weimer  of  Mission  Research  Inc. 
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Figure  1  A  Part  of  the  Model  Magnetic  Grid 
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Table  1:  Model  Prognostic  Variables 


Variable 

Name 

Units 

Description 

Comments 

Ion  densities 

particles' 

Local  (point)  volume  density  of  a  particular 

at  present  there  arc  7 

m3 

ion  species 

ions:  0\  H\  He\  N24 

02\  NO+,  N' 

Ion 

K.  degree 

Local  temperature  of  a  particular  ion 

same 

temperatures 

species 

Ion  velocities 

m/s 

Local  (point)  velocity  of  a  particular  ion 
species  along  the  magnetic  field  line 
passing  through  this  point 

same 

Electron 

temperature 

K.  degree 

Local  electron  temperature 

Electron 

m/s 

Local  (point)  velocity  of  electrons  along  the 

velocity 

magnetic  field  line  passing  through  this 
point 

Electron 

particles/ 

lx>cal  (point)  volume  density  of  electrons 

is  the  sum  of  all  local 

density 

m3 

Table  2:  Model  External  Variables 

ion  densities 

Variable 

Name 

Units 

Description  Comments 

Neutral 

particles/ 

Local  (point)  volume  density  -  at  present  there  are  7  neutrals: 

densities 

m 

of  a  particular  neutral  species  0,  02.  N2,  He.  H,  NO,  N. 

Neutral 

K. 

Local  temperature  of  all  same 

temperature 

degree 

neutral  species  (one  for  all) 

Neutral  zonal 

m/s 

Local  (point)  velocity  of  all  same 

velocity 

neutral  species  (one  for  all) 
in  the  zonal  direction  (east- 
west  eastward  is  positive) 

Neutral 

m/s 

Local  ( point)  velocity  of  all  same 

meridional 

neutral  species  (one  for  all) 

velocity 

in  the  meridional  direction 
( north-south,  northward  is 
positive) 
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Table  3:  Model  Diagnostic  Variables 


Variable  Name 

Units 

Description 

Comments 

ExB  zonal  velocity 

m/s 

Local  (point)  velocity  associated 
with  the  zonal  ExB  drift  of  the 
magnetic  field  line  passing  through 
this  point 

Needs  to  be  computed  from  empirical  ExB 
models 

ExB  meridional 
velocity 

m/s 

Local  (point)  velocity  associated 
with  the  meridional  ExB  drift  of  the 
magnetic  field  line  passing  through 
this  point 

same  as  abo\  e 

Photo  production 

Particles/s/m ' 

Number  of  particles  of  a  particular 
ion  species  produced  as  a  result  of 
photoiomzation  per  second  per  unit 
volume 

-  at  present  there  are  7  neutrals,  only  5  of  those 
can  be  photoiomzed  O.  O-  N?,  He.  N  Other 
ions  are  produced  via  chemical  reactions,  such 

as 

O*  +  H  —  H*  +  0. 

Chemical 

production 

Particles/s/m' 

Number  of  particles  of  a  particular 
ion  species  produced  as  a  result  of 
chemical  reactions  per  second  per 
unit  volume 

ITiere  are  2 1  chemical  reactions  at  the  present 
e.g .  0*  +  H  — ►  H1  +  O. 

Chemical  loss 

Particles/s/m’ 

Number  of  particles  of  a  particular 
ion  species  lost  as  a  result  of 
chemical  and  recombination 
reactions  per  second  per  unit  volume 

This  value  is  a  product  of  the  density 
(concentration)  of  the  ion  species  being 
destructed  and  the  chemical  loss  rate.  L 

Photoion  izatton 

rales 

l/S 

Coefficients  needed  to  compute 
photo  production 

-  at  present  there  are  7  neutrals  only  5  of  those 
can  be  photo-ionized  O.  02,  N>.  He.  N 

Chemical  reaction 

rates 

m’/s 

Coefficients  needed  to  oompute 
chemical  loss  due  to  electron 
exchange  reactions 

fhere  are  2 1  chemical  reactions  at  the  present 
e  g..  0+  +  H  —  1C  +  O 

Recombination 

reaction  rates 

l/s 

C  oefficients  needed  to  compute  loss 
due  to  recombination  chemical 

reactions 

There  arc  7  recombination  reactions,  e  g  .  (f  + 

e  —  O 

e  represents  an  electron 

Ion-neutral 

collision 

frequencies 

1/s 

Drag  on  a  particular  ion  particle  due 
to  collisions  with  a  neutral  species 

There  are  7  ions  and  7  neutrals,  therefore  it  is 
a  7x7  matr  ix  with  zero  diagonal 

Ion -ion  collision 
frequencies 

l/s 

Drag  on  a  particular  ion  particle  due 
to  collisions  with  a  different  ion 

species 

There  are  7  ions,  therefore  it  is  a  7x7  matrix 
with  zero  diagonal 

Ion  heating  rates 

Ion  thermal 

conductivities 

Electron  heating 

rates 

Electron  thermal 

conductivities 

J/m3/s 

J/K/m/s 

J/m3/s 

J/k/m/s 

Heating  due  to  Joule  heating, 
frictional  collisions  and  other 

processes. 

Heating  due  to  Joule  heating, 
frictional  collisions  and  other 

processes. 

Is  only  computed  for  three  major  ions,  04,  H\ 
He+ 

Is  only  computed  for  tliree  major  ions.  0\  H\ 
He* 
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2.1.1  Prognostic  Equations 

A  prognostic  equation  allows  one  to  estimate  a  particular  prognostic  variable  at  a  future  time.  The 
prognostic  variables  are  density,  velocity  and  temperature  for  ions  and  electron  density,  temperature 
and  velocity. 

These  equations  are  given  in  dipole  coordinates,  along  magnetic  flow  tubes.  Therefore,  there  is  only 
one  dependent  spatial  coordinate  corresponding  to  the  position  along  the  magnetic  flow  tube.  This  can 
be  a  non-dimensional  variable  q  or  a  dimensional  variables  =  q  •  Rc  (R^  is  the  radius  of  the  earth). 


2. 1 . 1 . 1  Continuity  Equation  For  Each  Ion  Species 

Numerical  solution  of  this  equation  should  generate  ion  density  N,  (/  +□/)  given  all  related  variables 
at  time  t. 


av, 

dt 


-b: 


ds 


+  N,-VV,+VN, 


■V, 


=  P-L.N. 


(I) 


where 


yV,  -  density  of  ion  i 

Vt  -  velocity  (aligned  with  the  magnetic  flow  tube)  of  ion  i 
s  =  q-Re 


b, 

P, 

L, 


=  yj\  +  3cos^(eccLat)  •  - 


V  ecc Radius , 

-  chemical  production  +  photochemical  production 

-  chemical  loss  rate 

•  N,  -  chemical  loss 


The  term 


vv± 


6  •  V[q  sin2 (eccLat)  •  (1  +  cos2  ( eccLat )) 
p  •  Rc  ■  (1  +  3  •  cos2  ( eccLat ))2 


(2) 


is  a  divergence  of  ExB  velocity  in  the  vertical  (and  meridional)  plane,  i.e.,  in  p  direction.  is  the 
value  of  ExB  meridional  drift  at  the  magnetic  equator  corresponding  to  a  particular  p. 

2. 1 . 1 .2  Momentum  Equation  For  Each  Ion  Species 

Numerical  solution  of  this  equation  should  generate  ion  velocity  Vf(t  +□/)  given  all  related  variables 
at  time  /. 
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I 


(3) 


V  = - 

/  A  _  Neutrals  A  Ions 

X  1  -  +  X  vij 

n  =  \  ./  =  I 


,  M, 

-g  sin  /  +  — — 


m 


i  11  i 


T,  dN,  |  T  |  a(r,  +  7-„) 
TV,  55  TV,.  55  55 


A'  Neutrals 


A  /««j 


+  X  cos  D  -Un  sin  D)cos  I  +  ^  ^;/F 


/=■• 


where 


TV,  -  density  of  ion  i 

Vi  -  velocity  (aligned  with  the  magnetic  flow  tube)  of  ion  i 


Vj  -  velocity  of  ion  j. 

s  =  q-  Re 


bs  =  y[l  + 3  cos2  (eccLat) 


R 


eccRadius 


mi  -mass  of  ion  i. 
k  -  Boltzmann’s  constant. 

Tt  -  temperature  of  ion  i. 

Te  -  electron  temperature 

Ne  -  electron  density. 

Un  -  zonal  neutral  velocity. 

Vn  -  meridional  neutral  velocity. 

Vm  -  ion-neutral  collision  frequency. 

Vtj  -  ion-ion  collision  frequency. 

g  -  acceleration  of  gravity. 

I  -  inclination  angle  for  this  flow  tube: 

.  _  2  cos  (eccLat) 

sin  I  =  .  ^  . 

yj\  +  3  cos 2 (eccLat) 
r  sin  (eccLat) 

COS  /  =  ^  - 

^1  +  3  cos2  (eccLat) 

D  -  declination  angle  for  this  flow  tube. 


2. 1 . 1 .3  Energy  Equation  for  Each  Ion  Species 

Numerical  solution  of  this  equation  generates  ion  temperature  Tf(t  +Dt)  given  all  related  variables  at 

time  t.  The  temperature  equation  is  solved  for  three  major  ions  (0\  Yt  and  He).  Temperatures  of  the 
remaining  ions  are  set  to  temperature  of  OX 


-kN, 

2 


+VyT  I  =  kNTti  — 


dt 


ds 


K 


d  (  dT,}  3 


dT 


-kN,Ti-'VV1+b* —  — -  +  —  kNtVibs  — L  +  Q  +  F  (4) 


ds  l.  ds 


ds 
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Where 


Q  and  F  -  are  different  types  of  heating  rates 
K  -  is  the  thermal  conductivity 


2.1. 1.3.1  Electron  Temperature  Equation 

It  is  similar  to  ion  temperature  equation,  except  that  the  conductivities  and  heating  rates  are  computed 
for  electrons. 


2.1 .1 .3.2  Electron  Density  Equation 

NumbcrOflons 

Ne  =  X  N<  (5) 

1*1 


2.1 .1 .3.3  Electron  Velocity  Equation 
Assumes  that  there  are  no  field-aligned  currents: 


A  'umberOflons 


I  kn, 


2.1.2  Fast  solution  of  tridiagonal  linear  systems 

After  discretization  of  the  prognostic  equations  we  obtain  a  system  of  linear  equations  which  has  the  tri¬ 
diagonal  structure.  A  square  matrix  A=(aii)  is  said  to  be  tridiagonal  if  Oy—0  for  all  pairs  (ij)  satisfying 
\i-j\  >  0.  Thus  in  the  rth  row,  only  au.i  and  a*  and  au+i  can  be  different  from  zero.  Three  vectors  can  be 
used  to  store  the  nonzero  elements. 

Using  standard  sweep  solver,  such  systems  can  be  solved  very  quickly.  In  fact,  the  solution  only 
requires  linear  time,  as  opposed  to  cubic  time,  in  the  general  case.  However,  if  the  matrix  is  ill- 
conditioned,  then  the  solution  can  be  inaccurate.  To  remedy  this,  we  explored  the  sophisticated 
pivoting  algorithms  of  LAPACK,  which  implement  careful  pivoting  and  are  robust  against  ill- 
conditioned  matrices  but  require  longer  execution  time.  We  determined  that  simple  sweep  solver  was 
sufficient.  The  code  also  implements  error  monitoring. 

2.13  Interpolation 

The  problem  of  interpolating  a  value  at  an  unknown  point  in  the  ionosphere  arises  in  two  important 
tasks  in  the  system.  The  first  problem  is  to  interpolate  information  about  the  neutral  atmosphere,  stored 
in  a  regular  geographic  grid  onto  each  point  in  the  magnetic  global  atmosphere  grid.  This  problem  is 
relatively  straightforward  because  of  the  regular  structure  of  the  grid  points. 

A  major  task  involves  interpolating  values  from  a  magnetic  (dipole)  grid.  This  is  required  for 
presenting  model  results  on  a  geographic  grid.  Also,  when  computing  slant  TEC,  one  needs  the  value 
along  an  almost  arbitrary  line  segment  passing  through  the  ionosphere  for  integration.  This  is  needed 
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for  assimilating  TEC  readings  from  ground  stations.  This  problem  is  more  difficult  because  there  is 
little  structure  of  the  density  of  the  points  to  exploit.  For  this  we  needed  more  sophisticated  geometric 
algorithms. 

We  describe  both  of  these  solutions  below. 

2.1. 3.1  Spherical  Interpolation 

We  use  “Spherical  Interpolation”  to  interpolate  a  value  from  a  regular  geographic  grid.  Because  this 
grid  locally  has  the  structure  of  a  Cartesian  grid,  it  is  straightforward  to  locate  the  cube  containing  the 
query  point.  The  interpolated  function  values  at  the  quety  point  can  be  obtained  by  using  trilinear 
interpolation. 

That  is,  suppose  we  have  a  point  at  position  (x,y,z)  located  inside  the  cube  in  Figure  2. 


<0.1,0  <1,1,0 


Figure  2:  Interpolation  from  a  Geographic  Grid. 

The  value  of  the  function  we  are  interpolating  at  position  (x,y,z)  within  the  cube  will  be  denoted  Vxyz 
and  is  given  by : 

Vxyz  =  Vooo(l  -x)(l  -  y)  ( 1  -z)  + 

VI00x(l-y)(l-z)  + 

V0,o(l-x)y(l-z)  + 

Vooi  (1  -x)(l  -y)z  + 

Vim  x(l  -y)z  + 

Von  (1  -x)yz  + 

Vi|0xy(l  -z)  + 

Vmxyz 


Although  technically,  the  “cubes”  in  the  geographic  grid  are  not  really  cubes  due  to  the  earth’s 
curvature,  the  density  of  the  grid  allowed  us  to  treat  them  as  cubes  with  negligible  error.  Indeed,  our 
algorithm  compared  well  with  the  interpolation  operations  in  MATLAB. 


2. 1 .3 .2  Interpolation  with  scattered  data 

When  interpolating  from  a  magnetic  grid,  the  problem  is  much  more  complicated,  because  there  is  not 
much  regularity  in  the  point  spacing. 

How  do  we  interpolate  a  function  when  we  know  the  function  at  points  whose  location  has  very  little 
structure  to  exploit?  We  make  use  of  a  fundamental  geometric  concept  known  as  a  triangulation  of  a 
set  of  points.  A  set  of  triangles  meeting  at  edges  and  vertices  whose  union  equals  the  domain  is  a 
working  definition  of  a  triangulation.  When  we  speak  of  triangulations  of  point  sets,  we  mean  that  we 
wish  to  form  triangulation  so  that  the  points  are  the  vertices  of  the  triangles,  with  the  boundary  edges 
being  the  edges  on  the  convex  hull  of  the  points.  (The  convex  hull  is  the  smallest  convex  set  containing 
the  point  set.) 

Having  a  triangulation  of  a  set  of  points  in  the  plane,  allows  us  to  locate  the  triangle  containing  the 
query  point,  and  then  writing  the  query  point  in  terms  of  the  barycentric  coordinates  of  that  triangle. 

We  can  then  use  the  barycentric  coordinates  to  weight  the  contribution  of  the  function  values  for  those 
three  points.  This  is  piece-wise  linear  interpolation. 

Because  a  set  of  points  in  the  plane  can  be  triangulated  many  ways,  what  is  the  best  triangulation?  This 
problem  has  been  considered  for  many  years  and  it  is  commonly  accepted  that  the  Delaunay 
triangulation  has  many  nice  properties  for  interpolation.  Further,  it  is  very  easy  to  compute.  In  our 
code,  we  used  the  implementation  given  in  Graphics  Gems  1.  Figure  3  shows  the  Delaunay 
triangulation  of  a  set  of  flow  tubes  of  nearly  equal  longitude,  projected  into  the  plane. 

While  two-dimensional  Delaunay  triangulations  are  well-suited  for  interpolating  two-dimensional 
functions,  the  three-dimensional  Delaunay  tetrahedralization,  is  not  as  easy  to  work  with.  The  biggest 
problem  is  that  they  are  more  difficult  to  compute.  Indeed,  they  would  consume  the  bulk  of  processing 
time  if  we  had  put  them  into  the  model.  To  retain  the  desirable  interpolation  properties,  we  exploited 
the  fact  that  a  set  of  flow  tubes  share  constant  magnetic  longitude.  So  to  interpolate  a  given  point  we 
found  the  two  slices  of  flow  tubes  surrounding  the  query  point.  We  then  projected  the  query  point  onto 
these  two  flow  tubes,  performed  barycentric  interpolation,  and  then  weighted  the  two  interpolations 
based  on  the  inverse  distance  from  the  query  point  to  the  interpolated  points  in  both  sets  of  flow  tubes. 
Thus,  rather  than  using  four  interpolating  points  in  the  magnetic  global  atmosphere  we  use  two  sets  of 
three  points  on  each  slice.  We  found,  that  by  interpolating  in  this  way,  we  were  nearly  within 
agreement  to  the  true  Delaunay  interpolation,  and  were  able  to  process  queries  much  faster. 
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Figure  3:  The  Delaunay  triangulation  of  a  portion  of  a  set  of  flow  tubes  with  the  same  magnetic  longitude. 


2.1.3.3  Line-of-Sight  Calculations 

To  compute  the  slant  TEC  from  the  model,  we  had  to  identify  the  coordinates  of  the  ground-station  and 
the  satellite  position.  We  then  discretized  the  line  segment  connecting  these  two  points,  and 
interpolated  the  electron  density  from  the  model  magnetic  grid  to  each  one  of  these  points  along  the 
segment.  Since  the  same  magnetic  grid  point  can  contribute  to  interpolated  values  at  several  segments, 
we  then  consolidate  all  weights  and  store  only  one  weight  per  grid  point.  These  weights  together  with 
indices  of  the  corresponding  points  constitute  observational  operator  in  the  Kalman  filter  section  of  the 
system. 

From  there,  we  used  the  composite-trapezoid  rule  to  integrate  the  density  from  each  of  these 
interpolated  values  and  obtain  a  model-simulated  slant  TEC.  We  explored  a  variety  of  step  sizes  when 
discretizing  the  line  segment,  and  chose  the  optimal  tradeoff  between  accuracy  and  running  time. 


2.1.4  Timing 

On  a  Pentium  Xeon  2.2GHz  computer,  1  model  time  step  requires  35  seconds  in  a  configuration  with 
106  grid  points  and  the  memory  footprint  is  1 .5  Gb. 


2.1.5  Tangent  Linear  and  Adjoint  Model 

We  consider  two  approaches  to  computing  evolution  of  background  error  covariance  in  the  system.  The 
first  is  described  in  Khattatov  et  al  (2000),  where  a  modified  tracer  mass  conservation  equation  is  used 
to  model  the  variance.  Off-diagonals  of  the  error  covariance  matrix  are  then  obtained  from  values  of  the 
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variance  using  an  assumption  of  Gaussian  de-correlation.  The  second  approach,  described  in  this 
section,  involves  obtaining  the  tangent-linear  of  the  underlying  model,  either  explicitly  as  a  sparse 
matrix,  or  implicitly  in  the  form  of  a  code  that  computes  a  product  of  the  linearization  matrix  and  an 
arbitrary  vector.  After  such  linearization  matrix  is  obtained  time  evolution  of  the  error  covariance  can 
be  explicitly  computed  as  described,  for  instance,  in  Khattatov  et  al  (2000).  After  some 
experimentation  we  decided  that  the  first  approach  is  more  practical,  while  the  second  should  be  used 
for  exploring  and  tuning  the  optimal  values  of  the  de-correlation  lengths.  We  therefore  describe 
calculations  of  the  linearization  matrix  with  respect  to  the  mass-conservation  equation  for  a  particular 
plasma  flow  tube. 

Let  us  assume  that  the  model  M  is  a  forward  time-dependent  discrete  propagator  that  accepts  the 
current  state  x,  (electron  or  ion  densities  throughout  the  model  domain  arranged  in  a  vector)  and  values 
of  several  atmospheric  drivers  p,  (e.g.,  level  of  solar  activity,  etc)  as  inputs  and  generates  state  estimates 
for  a  later  time: 


X,+A,  =M(p„  X,) 


(7) 


Both  quantities,  xtand  pt,  are  considered  to  be  model  parameters. 


Generally  speaking,  the  tangent  linear  of  the  model  M  is  simply  a  derivative  of  the  results  with  respect 
to  the  initial  conditions  or  input  parameters.  Note  that  M  is  a  non-linear  vector  function  and  therefore  its 
linearizations  (first  derivatives)  are  matrices: 


L 


X 


8M 

dx 


and  Lp 


8M 

dp 


(8) 


The  adjoint  of  the  propagator  model  is  simply  a  transposed  of  the  matrix  L  or,  a  way  to  compute  a 
product  of  the  transposed  and  an  arbitrary  vector.  The  linearization  matrix  describes  the  sensitivity  of 
the  model  with  respect  to  the  initial  conditions  and  the  adjoint  is  used  to  either  minimize  the  cost 
function  in  the  variational  assimilation  approach  or  to  solve  the  Kalman  filter  equations. 

These  linearization  matrices  can  be  obtained  in  several  ways: 

-  approximated  via  finite-differences  calculations,  i.e.,  by  introducing  small  changes  in  xt  and/or  pt  and 
computing  resulting  changes  in  x^f 


Ax,+a, 

Ax, 


(9) 


-  by  differentiating  the  actual  computer  code  (e.g.,  Fortran  or  C)  implementing  the  model  and 
generating  computer  code  for  direct  computation  of  L. 

-  by  analytical  differentiation  of  theoretical  equations  of  the  model  and  coding  the  results. 

The  first  approach  can  be  extremely  CPU  intensive  but  is  straightforward  to  implement.  The  other  two 
are  much  more  efficient  but  can  be  hard  to  implement  and  will  have  to  be  re-done  if  changes  are 
introduced  into  the  model. 

We  followed  the  second  approach  and  obtained  and  coded  explicit  calculations  of  both  these  matrices 
for  each  plasma  tube  in  the  model.  An  example  of  the  linearization  matrix  is  shown  in  Figure  4  for  0+. 
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Figure  4:  An  example  of  logarithm  of  absolute  value  of  linearization  matrix  for  one  plasma  tube. 

The  particular  plasma  tube  has  100  points.  Changes  in  ion  densities  at  the  next  time  step  can  be 
obtained  by  multiplying  this  matrix  by  a  vector  of  changes  in  ion  density  for  this  tube  at  a  previous  tie 
step.  The  asymmetries  in  the  matrix  reflect  the  fact  that  in  this  particular  case  the  prevailing  parallel 
transport  is  directed  from  the  top-left  end  of  the  tube  to  the  bottom-right. 

Once  the  linearization  matrix  is  obtained,  evolution  of  the  error-covariance  matrices  for  different  points 
on  the  same  tube  can  be  trivially  computed  as  follows: 

C(t  +  At)  =  LxC(t)Lrx 

2.2  Data  Assimilation 

The  Kalman  filter  is  a  mathematical  tool  traditionally  used  to  continuously  combine  model  predictions 
of  the  system  state  with  noisy  observations  in  order  to  come  up  with  a  less  uncertain  estimation  of  the 
system  state  and  error  covariances.  The  Kalman  filter  approach  adopted  here  resembles  that  of 
Khattatov  et  al  (2000).  Let  us  assume  that  model  estimates  of  electron  densities  at  all  grid  points  at  time 
t  are  arranged  in  a  vector  x  with  dimension  Nx.  Formally,  integration  of  the  model  M  can  be  written  as 

Xt+M  =  M0>  x/)  •  (10) 
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Let  vector  y  contain  observations  of  a  quantity  linearly  related  to  electron  densities  at  the  same  time.  In 
the  case  of  GPS  reference  station  data  such  quantities  are  slant  TEC  from  each  station  to  all  satellites  in 
view. 

The  connection  between  x  and  v  is  established  via  linear  interpolation  and  integration  observational 
operator  H  as  follows: 

y  =  H(x)  .  (11) 

Under  assumptions  of  linearity  and  Gaussian  statistics,  the  optimal  value  of  x  that  inverts  (1 1 )  given  a 
set  of  observations  y  and  model  estimates  of  x  is  given  by  the  Kalman  filter  [9]: 

x;  =x,  +K(y-Hx,)  (12) 

k  =  b,ht(Hb,ht+o  +  r)'1  .  (13) 

Here  B,  is  the  forecast  error  covariance  at  time  t.  O  is  the  error  covariance  matrix  of  the  observations 
and  R  is  the  representativeness  error  covariance  associated  with  errors  of  interpolation  and 
discretization.  Matrix  K  is  called  the  Kalman  gain  matrix. 

The  analysis  error  covariance  is  expressed  as: 

B?  =  B,  -B,Ht  (HB,Ht  +  0  +  R)  ‘HB,  .  (14) 

Once  inversion  of  (1 1 )  is  performed,  the  obtained  electron  densities,  x,“.  can  be  used  as  the  initial 
condition  for  the  model  M  to  predict  electron  densities  at  a  later  time  (beginning  of  the  next 
assimilation  window)  according  to: 

x,+A,  =  .  (15) 

Since  the  model  domain  contains  about  106  points,  direct  matrix  manipulations  described  by  (12),  (13), 
and  (14)  are  impossible  to  implement  even  with  modem  computing  capabilities.  As  is  customary  in 
large-scale  Kalman  filter  implementations  in  NWP  and  other  areas  of  atmospheric  sciences,  one  tracks 
only  the  evolution  of  the  diagonal  of  the  background  error  covariance  matrix  and  parameterizes  the  off- 
diagonal  elements.  If  we  also  assume  that  correlations  between  variations  of  electron  densities  at  two 
different  points  become  negligible  at  certain  distances  (i.e.,  assuming  compactly  supported  error 
covariance  models),  the  matrices  become  sparse  and  the  abovementioned  calculations  become  possible. 

We  argue  that  since  plasma  equations  are  solved  in  magnetic  coordinates,  in  order  for  the  background 
error  covariance  to  be  separable,  the  de-correlation  lengths  and  distances  between  points  need  to  be 
specified  in  magnetic  rather  than  geographic  coordinates  as  in  Hajj  et  al  (2003).  This  is  the  approach 
adopted  here. 

2.2.1  Kalman  Filter  and  Sparse-Background  Error  Covariance 

The  state  vector  X  consists  of  the  values  of  electron  density  Ne  in  every  node  of  the  computational 

grid.  Given  the  size  of  the  state  vector  X  (about  a  million  entries),  the  model  covariance  matrix  B  (up  to 
million  by  million  entries),  and  the  observational  matrix  H  (up  to  600  by  a  million  entries),  the  exact 
implementation  of  the  abovementioned  equations  is  out  of  the  question. 
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The  model  covariance  matrix  B  contains  the  information  about  the  mutual  correlation  of  the  values  of 
electronic  densities  in  every  grid  cell  with  each  other.  Consequently,  its  size  is  equal  to  the  number  of 
grid  points  squared.  Given  the  large  size  of  this  matrix,  it  is  not  feasible  to  properly  update  all  of  its  off- 
diagonal  elements  during  every  Kalman  Filter  step.  Instead,  only  the  evolution  of  the  variance  (i.e. 
diagonal  of  the  covariance  matrix  B)  is  tracked  during  the  computation,  and  the  off-diagonal  elements 
are  parameterized  by  the  diagonal  elements  of  B  and  are  simply  rescaled  every  time  according  to  the 
recomputed  variance.  The  time  update  is  also  performed  for  the  variance  only. 

The  off-diagonal  elements  Bt  correspond  to  the  correlations  between  the  corresponding  nodes  of  the 
magnetic  grid.  I^t  /  and  j  correspond  to  the  two  nodes  of  the  magnetic  grid  with  the  magnetic 
coordinates  loni ,pi^qi  and  Ion }  ,prqr  respectively.  The  correlation  Cort  ;  between  these  two 
nodes  is  determined  as  follows: 


Cor,j  =  exP(- 


-)  *exp( 


-)  *  exp( — 


Ion 


and  the  corresponding  entry  in  the  covariance  matrix  B  is  determined  as 

B,j=^^Coru- 

To  make  the  amount  of  nonzero  B  elements  manageable,  all  correlations  Cor(  that  are  less  then  the 

threshold  of  0.01,  are  discarded.  The  correlations  between  a  particular  fixed  grid  point  and  all  other 
grid  points  with  the  same  value  of  magnetic  longitude  and  with  the  same  value  of  p  are  presented  in 

Figure  5. 
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Figure  5:  An  example  of  q-p  and  magnetic  longitude  -  q  cross-sections  of  background  error  covariance 
between  a  fixed  point  in  the  model  grid  and  all  other  grid  points. 
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The  variance  Bt  t  is  recomputed  after  every  Kalman  Filter  analysis  step  as  follows: 


Ba  =  (/  -  KH )  Bh ,  where  Ba  is  the  resulting  variance,  Bh  is  the  variance  before  the  Kalman  Filter 
analysis  step,  I  is  the  unit  matrix  and  K  is  the  Kalman  gain  matrix,  defined  as 
K  =  BH '  I HBH 7  +  O  +  R )  ,  with  O  being  a  covariance  matrix  of  the  observational  errors,  and  R 
being  a  regularizer  matrix. 

The  time  update  of  the  covariance  matrix  B,  Bl+Al  =  F  (  Ba  )  is  also  performed  only  on  the  variance, 

with  the  off-diagonal  elements  of  B  rescaled  accordingly  afterwards,  as  shown  above.  The  variance 
update  is  performed  as  follows: 


B, 


t+ to 


=  By+Ql+J-'=B;-+a'Ne; 


Here  a  is  a  relative  model  error  accumulated  between  successive  analyses  and  Ne,  is  an  electron 
density  in  the  grid  point  i. 

Special  care  needs  to  be  taken  regarding  the  speed  and  memory  requirements  of  the  implementation  of 
the  Kalman  Filter  equations.  A  set  of  custom  data  structures  has  been  designed  to  store  the  sparse 
observational  matrix  H,  various  lookup  tables  for  the  correlations  between  model  grid  points,  and 
various  intermediate  matrices.  In  theory,  it  is  possible  to  precompute  all  the  correlations  during  the 
initialization  phase.  However,  given  that  each  correlation  requires  up  to  1 0  bytes  of  storage  (at  least  4 
bytes  for  the  correlation  value  and  up  to  6  bytes  to  store  information  about  what  grid  points  it 
corresponds  to),  the  amount  of  memory  required  to  store  correlations  would  be  prohibitively  large.  In 
addition,  the  experiments  showed  that  the  time  to  retrieve  the  precomputed  correlations  can  be 
significant,  given  the  large  size  of  the  lookup  table  and  the  totally  random  manner  in  which  correlations 
need  to  be  accessed.  Consequently,  the  correlations  are  mainly  computed  on  the  fly. 

Such  an  approach  could  also  be  prohibitively  expensive  if  performed  in  a  straightforward  fashion.  Of 
course,  not  all  elements  of  B  are  computed  each  time,  but  only  those  that  have  corresponding  nonzero 
entries  in  the  observational  matrix.  Physically  it  means  that,  during  each  assimilation  step,  we  only  need 
to  compute  the  correlations  between  the  model  computational  cells  intersected  by  line  of  sights  between 
receivers  and  satellites,  and  all  other  computational  cells.  In  addition  to  the  abovementioned 
simplification,  the  algorithm  to  compute  the  product  of  the  covariance  matrix  B  and  the  observational 
matrix  H  takes  advantage  of  the  properties  of  the  magnetic  grid  and  the  “fractal”  structure  of  the 
covariance  matrix  B.  The  algorithm  also  takes  advantage  of  the  fact  that  correlations  are  considered 
only  if  their  value  is  less  than  a  certain  threshold,  and  uses  an  intelligent  way  to  make  a  quick  estimate 
of  a  correlation  without  determining  its  actual  value.  To  speed  up  the  computation,  several  lookup 
tables  are  created  at  start-up  that  contain,  in  a  concise  form,  the  information  about  what  ranges  of 
elements  of  the  covariance  matrix  B  are  nonzero.  This  information  is  used  in  conjunction  with  the 
bounding  box  for  a  particular  observation  (i.e.  a  line  of  sight  between  a  receiver  and  a  satellite)  to 
quickly  determine  the  (hopefully)  small  subset  of  the  correlations  that  actually  need  to  be  computed. 

2.2.2  Slant  TEC  determination  from  the  IGS  Global  Positioning  System  data. 

The  IGS  (International  GPS  Service)  Global  Positioning  System  data  is  utilized  to  compute  the  receiver 
and  satellite  coordinates,  and  to  determine  the  total  slant  TEC  value  along  the  lines  of  sight  between 
receivers  and  satellites.  The  algorithm  to  determine  the  latter  is  similar  to  that  described  in  Blewitt 
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(1990).  The  dual-frequency  GPS  pseudorange  data  from  the  GPS  ground  stations  include  carrier  phase 
and  code  measurements,  both  of  which  have  to  be  utilized  in  determining  slant  TEC  values.  The  code 
data,  which  accuracy  is  in  the  meter  range,  is  used  to  determine  the  initial  integer  number  of  phase 
cycles  between  a  satellite  and  a  ground  station.  The  phase  data  is  much  more  accurate,  with  accuracy 
in  the  millimeter  range,  but  the  phase  data  ambiguity  regarding  the  integer  number  of  cycles  between 
satellites  and  receivers  has  to  be  accounted  for.  This  number  is  initially  determined  with  the  help  of  the 
code  data,  and  it  remains  constant  for  as  long  as  the  receiver  maintains  a  lock  on  the  signal.  However, 
the  so  called  “cycle  slips,”  or  sharp  discontinuities  in  the  carrier  phase  data,  that  happen  when  a  receiver 
fails  to  maintain  a  lock  on  the  satellite  signal,  have  to  be  accounted  for.  This  is  accomplished  by 
maintaining  a  cache  of  approximately  40  minutes  worth  of  previous  data. 

2.23  Bias  Determination 

Dual-frequency  GPS  receivers  normally  have  large  differential  code  biases  (DCB).  Small  differences  in 
the  LI  and  L2  hardware  channels  of  dual-frequency  receivers  result  in  additional  differential  signal 
delays.  If  ignored,  these  biases  will  skew  the  estimated  slant  TEC  measurements. 

Quite  often  these  biases  are  as  large  as,  or  even  larger  than,  the  useful  signal.  Typical  measured  TECs 
are  in  the  range  of  0-100- 1016  electrons  per  m2  (or  1-100  TEC  units,  TECU).  Receiver  biases  expressed 
in  TEC  units  that  we  have  observed  can  easily  reach  ~  ±50  TECU.  These  biases  can  also  change  (albeit 
slowly)  and  are  the  effect  of  aging  equipment  or  varying  environmental  conditions.  Clearly,  for  the 
reference  station  data  to  be  useful  such  biases  need  to  be  estimated,  monitored,  and  removed  from  the 
measurements. 

A  simple  way  to  estimate  bias  is  to  assume  zero  ionospheric  electron  content  at  night.  The  ionosphere 
however,  can  easily  vary  by  ~  10  TECU  through  the  night  and  this  introduces  no  insignificant  errors  in 
the  estimated  bias. 

A  more  sophisticated  way  relies  once  again  on  the  “thin  shell”  approximation.  In  a  “nutshell,”  vertical 
TEC  is  simultaneously  estimated  at  a  number  of  ionospheric  piercing  points  between  station  i  and 
satellite  j.  Then  it  is  assumed  that  the  TEC  variations  between  two  geographic  locations  are  correlated. 
Typically,  a  Gaussian  correlation  function  is  used  with  de-correlation  lengths  of  a  few  hundred 
kilometers  specified  in  geographic  coordinates.  The  sum  of  the  unknown  bias  and  the  unknown  ‘true’ 
slant  TEC  constitutes  the  measured  slant  TEC[y].  Unknown  biases  and  ‘true’  TECs  can  be  estimated 
using  linear  Kalman  filter  starting  from  a  reasonable  background  guess.  After  enough  data  is 
accumulated  and  if  the  error  covariance  model  and  background  guess  are  correct,  uncertainties  in  both 
unknowns  will  decrease.  A  simplification  of  this  approach  assumes  that  nearby  true  TECs  are  the  same 
if  they  are  located  within  a  certain  distance.  This  way,  if  enough  data  are  available,  one  can  write  a 
system  of  linear  equations  relating  measurements  and  unknowns.  If  biases  from  one  or  more  stations 
are  known  a  priori  the  system  can  be  solved  and  biases  estimated. 

The  quality  of  bias  estimation  clearly  depends  on  the  validity  of  the  thin  shell  approximation  and,  in 
turn,  directly  affects  the  quality  of  the  differential  ionospheric  corrections,  which  are  as  well  influenced 
by  the  errors  of  the  “thin  shell”  approximation  and  errors  of  the  interpolation. 

Formally,  slant  TEC  measured  by  a  particular  station  is  a  sum  of  ‘true’  and  unknown  slant  TEC  and 
receiver  and  satellite  biases: 

=  y,rue  +br+b,  (16) 
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Depending  on  the  local  time,  the  magnitude  of  these  biases  is  often  significantly  larger  than  the  actual 
slant  TEC.  In  order  for  the  data  to  be  relied  upon  these  biases  need  to  be  determined  and  accounted  for. 
This  can  be  accomplished  by  augmenting  vector  x  with  both  types  of  biases.  Since  the  size  of  x  is  ~  1 06, 
the  number  of  satellite  biases  is  ~  30,  and  the  number  of  IGS  stations  we  currently  use  is  ~  100,  this 
does  not  lead  to  any  significant  increase  of  the  number  of  unknowns. 

In  principle,  both  satellite  biases  and  station  biases  can  be  continuously  computed  using  this  method.  At 
this  stage  in  the  development  process  we  chose  to  use  fixed  broadcast  satellite  biases  and  only 
determine  receiver  DCBs.  This  can  be  justified  by  noting  that  satellite  biases  are  perceived  to  be 
generally  better  known  than  the  receiver  biases  and  they  vary  at  slower  rate.  We  do  plan  to  add 
continuous  satellite  bias  calculations  in  the  future. 

Figure  6  illustrates  bias  adjustment  results  for  two  different  stations  along  with  DCBs  determined 
independently  by  two  IGS  centers. 


Figure  6:  Examples  of  bias  estimation  obtained  in  the  system  (black  line)  are  shown  next  for  two  different  stations  together  with  bias 
estimation  results  from  other  institutions  when  available  (JPL  -  magenta;  CODE  -  green). 


A  sample  of  the  derived  biases  is  shown  in  Table  4  for  several  IGS  reference  stations  along  with  IGS 
Ionosphere  Working  Group  computed  biases.  As  one  can  see  in  most  cases  (and  this  is  a  typical  for  the 
— 145  IGS  stations  we  used  so  far)  the  two  sets  of  biases  agree  with  each  other  to  within  several  TEC 
units.  Several  cases  with  large  differences  are  still  present,  however.  The  ultimate  validation  of  these 
biases  will  come  from  systematic  comparisons  of  the  assimilation  analysis  produced  in  the  system  with 
independent  data. 
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Table  4:  A  Sample  of  Estimated  Receiver  Biases 


Station  Name 

Estimated  biases,  TEC  units 

IGS  Ionosphere  Working  Group 
biases  from  JPL,  TEC  units 

ALBH 

-23 

-28 

BAHR 

29 

29 

CAS1 

-18 

-18 

DAVI 

-45 

-38 

FAIR 

22 

18 

GUAM 

42 

40 

HOFN 

10 

13 

use 

35 

25 

JPLM 

-7 

-15 

KARR 

2 

0.6 

MADR 

10 

17 

NLIB 

-8 

-14 

OHI2 

7 

4 

PERT 

19 

18 

It  normally  takes  several  diurnal  cycles  for  the  bias  calculations  to  stabilize.  After  about  5-10  days 
oscillations  in  the  biases  decrease  significantly.  The  described  bias  determination  approach  appeared  to 
work  reasonably  well  for  a  period  of  time.  However  sometimes  obtained  biases  show  large  diurnal 
variations  on  the  order  of  10  TEC  units.  We  attribute  this  to  the  model  having  a  time-dependent  bias  as 
well.  We  are  currently  developing  s  method  for  simultaneously  estimating  model  biases  together  with 
data  biases  in  the  same  Kalman  filter.  Until  this  work  is  completed  we  recommend  using  externally- 
derived  DCBs  in  operational  implementation. 

2.3  Data 

There  are  two  types  of  input  data  necessary  for  system’s  operations:  static  that  does  not  change  with 
time  and  dynamic  that  does  change  with  time  and  therefore  needs  to  be  constantly  updated. 

23.1  Static  Input  Data 

Static  data  includes  the  following: 

•  Solar  EUV  flux  at  different  wavelengths  that  is  later  scaled  in  the  model  by  the  value  of  the 
FI 0.7  flux. 

•  Absorption  cross-sections  for  several  neutral  gases. 

•  Photoionization  cross-sections. 
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•  Coefficients  of  the  World  Magnetic  Model  2000  (WMM  2000; 

http://www.ngdc.noaa.gov/segAVMM/DoDWMM.shtinn. 

•  Pre-computed  monthly  3-D  fields  of  neutral  composition  and  temperature  from  MS1SE  model. 

(http://nssdc.gsfc.nasa.gov/space/model/atmos/msise.html ) 

•  Pre-computed  monthly  3-D  fields  of  neutral  winds  from  HWM  model  (http://uap- 
www.nrl.navy.mil/models  web/hwm/hwm  home.htm). 


232  Dynamic  Input  Data 

2.3.2. 1  Solar  Fluxes 

One  of  the  most  important  drivers  for  the  described  system  is  the  current  level  of  solar  activity.  These 
data  (e.g.,  fl0.7  solar  flux)  are  obtained  every  1 5  minutes  form  the  NOAA  Space  Environment  Center, 

http://www.sec.noaa.gov. 

2.3 .2.2  GPS  Reference  Station  Data 

Navigational  RINEX  files  that  contain  satellite  orbit  information  are  obtained  from  IGS  network 
stations  archived  at  NASA’s  Crustal  Dynamics  Data  Information  System  (CDDIS)  with  either  hourly 
or  daily  delays  (http://cddisa.gsfc.nasa.gov/gps  datasum.html).  We  currently  use  1  -day  delayed  data  to 
simplify  operations. 

Daily  or  near-real-time  reference  station  observational  RINEX  data  files  are  also  obtained  from  the 
CDDIS  archive  as  compact  RINEX  and  uncompressed  for  processing  as  needed 

We  obtain  IGS  Ionosphere  Working  Group  produced  satellite  and  receiver  biases  in  IONEX  files  from 
ftp://cddisa.gsfc.nasa.gov/pub/gps/products/ionex.  These  files  are  only  available  with  about  a  7  day 
delay. 

2.3.2.3  Data  Handling  Scripts 

In  order  to  handle  updating  dynamic  data  files  we  developed  several  redundant  data-fetching  scripts 
residing  on  Linux  server  that  continuously  bring  in  (via  ftp  or  http)  and  archive  the  data.  Data  latency  is 
monitored  and  recorded  and  operator  can  be  alerted  to  possible  problems. 

2.4  Data  Processing  and  System  Operations 

An  overview  of  data  processing  is  given  in  Figure  7.  Data  collection  is  independent  from  the  model  and 
assimilation  system  operation,  if  data  is  delayed  for  whatever  reason,  a  warning  is  generated  and  the 
system  operates  with  the  most  recent  solar  data  and  no  data  assimilation. 

The  model  time  step  is  1  minute.  Model  time  is  synchronized  with  real-time.  So  occasionally  the  model 
has  to  stay  idle  to  maintain  synchronization.  At  the  beginning  of  each  time  step,  the  most  recent  solar 
activity  data  is  obtained  and  used  to  project  current  3-D  distributions  of  electron  densities  one  step  into 
the  future.  The  system  then  estimates  locations  of  all  satellites  and  computes  the  slant  TECs  from 
several  user-set  locations  on  Earth  to  all  satellites  in  view.  Results  are  displayed  on  the  web  page. 
Additionally,  a  profile  of  electron  densities  corresponding  to  the  location  of  an  operational  dynasonde 
at  Utah  State  University’s  Bear  Lake  Observatory  is  extracted  and  archived  for  comparisons  with  the 
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ionosonde  measurements.  NOAA’s  group  accesses  the  profile  and  displays  real-time  comparisons  on 
their  web-page.  Every  hour  the  model  3-D  fields  of  electron  and  ion  concentrations,  densities,  and 
velocities  are  saved  on  disk  and  archived  for  up  to  five  days. 

When  the  system  operates  in  the  assimilation  mode,  after  model  calculations  are  completed  at  each 
time  step,  phase  and  code  measurements  from  observational  RINEX  files  for  each  receiver-satellite  pair 
are  read  in  and  combined  to  obtain  a  vector  of  smoothed  measured  slant  TEC.  This  vector  corresponds 
to  vector  y  in  the  Kalman  filter  equations.  The  resulting  analyzed  electron  density  field  is  partitioned 
between  different  ions  in  the  model  and  these  become  initial  conditions  for  the  next  model  time  step. 


Figure  7:  Data  Flow  in  the  System. 
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2.5  Ensemble  Forecasting  and  Targeted  Observations 

We  explored  the  strategy  for  finding  locations  of  “most  important”  observations  of  the  electron  density 
in  the  ionosphere  using  the  ensemble  prediction  simulations. 

The  concept  of  adaptive  observations  aims  to  find  specific  geographical  areas  and  time  periods  that 
need  to  be  observed  for  obtaining  the  best  future  forecasts.  The  ultimate  goal  of  adaptive  (or  targeted) 
observations  is  to  decrease  the  forecast  errors  by  placing  sensors  (or  planning  special  observational 
modes)  in  regions  where  additional  data  are  expected  to  improve  the  forecast.  In  the  atmosphere  and 
ocean  studies,  the  strategy  of  adaptive  observations  is  based  on  the  recent  progress  in  the  understanding 
and  analysis  of  predictability  and  sensitivity  of  geo-systems  that  were  achieved  through  extensive 
applications  of  the  ensemble  based  techniques.  In  order  to  identify  the  most  sensitive  regions  for 
adaptive  observations,  information  about  the  time  evolution  of  uncertainty  in  the  initial  state,  boundary 
conditions  and  the  state  drivers  (forcing  terms)  should  be  analyzed  {Farrell  and  Ioannou ,  2004). 

The  Initial  State  (IS)  uncertainty-based  targeting  exploits  time  evolution  of  the  analysis  and  forecast 
errors  to  identify  the  most  important  regions  for  adaptive  observations.  Several  objective  techniques 
have  been  developed  which  promise  to  predict  the  optimal  sites  for  additional  observations.  The 
differences  between  forecasts  based  on  the  routine  plus  adaptive  observations  and  only  routine  data 
provide  evaluation  of  the  impact  of  adaptive  sensors.  The  breeding  method  of  Kalnay  and  Toth  (1997), 
the  Singular  Vector  (SV)  techniques  of  Buzzia  etal.{\  998),  and  ensemble  of  the  replicated 
(‘perturbed’)  observations  {Lorenz  and  Emanuel,  1998)  have  been  compared  and  discussed  in  a  number 
of  previous  studies.  Each  of  these  techniques  has  certain  attractive  features.  For  instance,  the  breeding 
method  can  provide  information  on  the  unstable  regions  of  the  finite  amplitude  perturbation  without 
additional  software  development  as  needed  for  application  of  SV  techniques.  Lorenz  and  Emanuel 
(1998)  concluded  that  the  ensemble  based  techniques  with  the  perturbed  observations  provide  more 
consistent  strategy  for  the  adaptive  targeting  in  the  simple  scalar  forecast  model.  In  the  Naval  Research 
Laboratory  (NRL)  the  adaptive  observation  approach  is  based  on  singular  vector  and  sensitivity 
functional  algorithms  obtained  with  the  aid  of  the  adjoint  versions  of  global  and  regional  operational 
models  {Baker  and  Daley.  2000). 

For  the  space  weather  studies  the  ensemble  based  data  analysis  techniques  and  adaptive  targeting  are 
relatively  new  disciplines.  The  direct  application  of  the  results  from  the  numerical  weather  prediction 
studies  might  be  erroneous  because  of  different  temporal  and  spatial  scales  of  physics  and  chemistry. 
For  example,  the  most  energetic  ionospheric  waves  present  ageostrophic  type  of  oscillations.  The 
strong  global  diurnal  variations  create  the  major  impact  on  the  redistribution  of  the  neutrals,  ions,  and 
electrons  in  the  upper  atmosphere.  In  comparison  with  the  tropospheric  techniques  that  mostly  intend  to 
target  positions  of  unstable  zones,  in  the  ionosphere  we  need  also  to  link  position  and  timing  of  the  fast 
unstable  dynamics.  The  observed  sunset  equatorial  scintillations  present  a  good  example  of  the  ‘time- 
space’  targeting  of  the  unstable  regions  in  the  ionosphere.  Below  we  describe  our  approach  and  results 
for  adaptive  observations  of  the  electron  density  using  the  ensemble-based  framework  with  the 
uncoupled  ionosphere  plasma  model. 

The  uncoupled  feature  of  the  model  formulation  assumes  the  specification  of  the  large  set  of  the  input 
parameters  including  the  neutral  winds  and  compositions,  ExB  drift  velocities  from  the 
electrodynamics,  etc.  The  approximate  (climatological/empirical)  knowledge  of  these  input  model 
drivers  can  be  considered  as  an  important  source  of  the  forecast  errors  of  electron  and  ion  densities. 

In  the  data  assimilation  the  quality  of  the  initial  conditions  contribute  to  determining  the  quality  of  the 
forecast  state.  The  traditional  formulation  of  the  adaptive  observation  is  to  study  the  model  sensitivity  to 
the  perturbation  of  initial  states  and  find  the  regions  where  the  forecast  errors  grow  quickly  and  linear 
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perturbation  of  the  control  forecast  state  may  be  unstable.  In  other  words,  by  placing  the  additional 
sensors  in  those  unstable  zones  we  hope  to  collect  the  most  ‘desirable’  information  that  can  help  to 
improve  the  quality  of  the  consecutive  analysis.  A  breeding  vectors  approach  is  used  to  find  locations 
of  dominant  perturbation  modes  using  the  ensemble  forecast  simulations.  For  the  ionosphere  plasma 
model  the  bred  vectors  are  produced  as  follows:  a  control  (base)  run  B  is  started  from  suitable  initial 
conditions,  then  the  perturbation  run  P  is  begun  from  the  same  initial  state  but  with  added  random 
(normally  distributed)  perturbations  of  the  electron  density.  Perturbed  forecast  P  proceeds  for  1  hour 
then  the  perturbations  defined  as  the  differences  between  P  and  B  at  all  model  grid  points  are 
renormalized  using  a  prescribed  norm.  In  case  of  the  uncoupled  plasma  forecast  model  we  use  the 
initial  root  mean  squared  variance  of  the  stochastic  perturbations  of  the  electron  density.  This 
renormalization  is  repeated  every  hour  until  the  perturbations  have  formed  the  relatively  well-defined 
perturbed  modes.  Several  perturbed  simulations  form  our  ensemble  prediction  system.  In  this  report  we 
use  various  sizes  of  ensembles  (5,  10,  and  20)  with  variable  strength  of  initial  perturbations  (5%,  10%, 
and  20%).  The  size  of  the  ensemble  determines  the  local  dimensions  of  the  bred  vectors  V.  In  the 
breeding  approach,  renormalization  of  perturbations  is  necessary  to  avoid  saturation  of  resolved 
unstable  oscillations.  It  is  worthy  to  highlight  that  the  saturation  of  the  noisy  oscillations  irrelevant  to 
the  forecast,  is  naturally  controlled  by  the  nonlinear  breeding. 

The  covariance  matrix  C  —  VV  can  be  considered  as  a  proxy  for  the  local  covariance  ensemble 
covariance  matrix.  Patil  et  al.  (2001)  showed  that  the  C-matrix  can  be  used  to  evaluate  the  ‘actual’  local 
finite-time  Bred  Vector  (BV)  dimension  employing  the  distribution  of  the  singular  vectors  of  C-matrix. 
In  the  numerical  weather  prediction  the  low  BV-dimensions  as  well  as  a  structure  of  the  singular 
vectors  of  C-matrix  can  indicate  the  most  unstable  subspace.  As  pointed  out  in  the  previous  studies 
(e.g.,  Kalnay  2000)  the  bred  vectors  by  construction  are  closely  related  to  Lyapunov  Vectors  (LV)  and 
Singular  Vectors  (SV)  of  the  linear  tangent  model.  The  bred  vectors  share  some  properties  of  the 
leading  LVs  including  the  following:  a)  BVs  are  independent  of  the  norm  used  to  define  the  size  of  the 
perturbation;  b)  BVs  are  independent  on  the  length  of  the  rescaling  time  window  as  long  as  strength  of 
the  perturbation  remains  linear.  The  most  significant  differences  between  BVs  and  LVs  can  be  invoked 
from  the  BV  definition  or  construction.  BV  are  finite-amplitude,  finite-time,  and  because  they  are  not 
globally  orthogonalized,  they  bear  the  local  properties  in  space.  In  the  next  section  we  present  our 
results  on  the  predictability  and  strategy  of  adaptive  targeting  for  electron  density  observations  and 
forecast  using  the  described  above  the  ensemble  breeding  technique. 


3.  RESULTS  AND  DISCUSSION 

3.1  Model  Results 

In  this  section  we  give  a  brief  overview  of  the  model  results.  Since  these  results  do  not  have  much 
scientific  value,  only  a  brief  description  is  provided  to  demonstrate  that  the  model  produces  physically 
realistic  output. 

Figure  8  and  Figure  9  show  typical  vertical  cross-sections  of  densities  (left),  parallel  velocities  (middle) 
and  temperatures  for  electrons  and  all  model  ions  at  a  fixed  magnetic  longitude. 

Typical  distributions  of  ionospheric  electron  densities  obtained  in  the  model  are  shown  in  Figure  10  as 
a  three-dimensional  isosurface  superimposed  over  the  Earth.  Concentrations  of  electrons  inside  this 
isosurface  are  larger  than  51012  electrons/m2;  electron  densities  outside  are  less  than  that.  The  radial 
direction  has  been  exaggerated  to  show  detail. 


25 


The  distribution  clearly  shows  equatorial  depletion  and  low-latitude  enhancements  that  arise  due  to 
vertical  transport  of  ions  (across  magnetic  field  lines)  at  the  equator  followed  by  descent  along  the  field 
lines  due  to  gravity.  This  is  a  well-known  feature  and  the  model  realistically  reproduces  it. 

Figure  1 1  presents  latitude-longitude  distribution  of  a  typical  vertical  TEC  in  the  model.  Note  that  since 
high-latitude  ionization  processes  are  currently  turned  off,  this  plot  does  not  show  enhanced  electron 
densities  in  the  Polar  regions. 
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Figure  8:  Examples  of  instantaneous  model  prognostic  variables  for  electrons  and  major  ions  shown  at  a  fixed  magnetic  longitude 

Real  time  system  data  are  being  currently  downloaded  to  NOAA  and  compared  operationally  with  the 
Bear  Lake  observatory  dynasonde  data.  Figure  12  shows  48  hours  of  comparisons  between  model¬ 
generated  TEC  at  the  dynasonde  location  and  experimental  data.  Detailed  information  on  various 
characteristics  of  Bear  Lake  ionosonde  can  be  accessed  here: 

http://www.ngdc.noaa.gov/stp/IONO/Dvnasonde/. 
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General  features  of  the  vertical  TEC  evolution  are  similar  in  both  time  series.  Moreover,  for  extended 
periods  of  time  the  difference  between  the  observations  and  the  model  is  only  about  a  few  TEC  units, 
which  is  rather  encouraging.  However,  the  model  clearly  does  not  show  as  rapid  decay  of  electron 
densities  at  sunset  as  the  ionosonde  data  suggest.  One  possible  reason  for  this  discrepancy  is  an  absence 
of  one  or  more  chemical  recombination  mechanism  in  the  model  responsible  for  fast  conversion  of  ions 
to  neutral  particles.  Another  possible  reason  is  that  dynasonde  measurements  extend  only  to  the  F-layer 
maxima  (about  300  km).  A  semi-empirical  correction  needs  to  be  added  to  the  derived  vertical  TEC  to 
account  for  missing  part  of  the  profile  (for,  instance,  via  vertical  extrapolation).  More  research  is 
needed  to  properly  explain  this  discrepancy.  Also,  the  model  nighttime  TEC  is  generally  lower  than 
dynasonde  measurements.  This  is  likely  to  be  due  to  missing  nighttime  ionization  source.  Realistic 
nighttime  ionization  will  be  included  in  the  model  in  the  near  future. 
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Figure  9:  Examples  of  instantaneous  model  prognostic  variables  for  minor  ions  shown  at  a  fixed  magnetic  longitude. 
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Figure  10:  Examples  of  model -generated  3-D  instantaneous  iso-surfaces  of  electron  density. 
The  radial  direction  is  exaggerated  to  show  detail. 
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Figure  1 1 :  Examples  of  model-generated  total  vertical  electron  content. 
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Figure  12:  Time  series  of  Bear  Lake  dynasonde  measurements  (red  dots)  and  model  simulated  vertical  TEC  (solid  blues). 
Courtesy  of  F.T.Berkey  (USU's  Bear  Lake  Observatory),  J.W.Wright  and  N.A.Zabotin  (both  at  CIRES/NGDC). 


Figure  13:  Electron  density  profile  measured  by  Bear  Lake  dynasonde  (red)  and  simulated  in  the  model  (blue). 
Courtesy  of  F.T.Berkey  (USU’s  Bear  Lake  Observatory),  J.W.Wright  and  N.A.Zabotin  (both  at  CIRES/NGDC). 


Figure  13  shows  a  vertical  profile  of  electron  density  measured  by  the  dynasonde  and  simulated  by  the 
model.  Once  again,  while  generally  the  plots  showed  the  same  features,  some  quantitative  discrepancy 
remains. 

The  polar  regions  are  characterized  by  a  variety  of  distinct  physical  processes  due  to  solar  wind  and 
interplanetary  magnetic  field  penetrating  the  magnetosphere.  Due  to  the  extra  complexity  involved  in 
modeling  these  processes  we  decided  to  leave  some  of  these  processes  out  and  instead  concentrate  on 
developing  the  system  at  low  and  middle-latitudes.  For  example,  auroral  precipitation  and  heating 
processes  are  not  included  in  the  model.  Nevertheless,  due  to  importance  of  the  high-latitude  regions 
we  included  a  description  of  the  most  fundamental  phenomena  with  the  goal  to  fill  in  the  missing 
processes  later. 

The  model  grid  extends  to  about  85  degrees  for  both  Northern  and  Southern  Hemispheres.  Earth’s 
magnetic  field  is  computed  in  the  polar  regions  (Figure  14)  and  Weimer  high-latitude  ExB  drift  model 
is  used  to  obtain  electric  potential  (Figure  15)  and  then  ExB  drift  vector  (Figure  16).  Three  components 
of  this  vector  are  later  used  for  computing  particle  transport  using  a  semi-Lagragian  approach.  An 
example  of  trajectories  for  three  particles  is  shown  in  Figure  17. 


29 


Figure  15:  The  x-  and  y-component  of  the  Weimer  electric  potential  at  the  North  Pole 


Figure  16:  Calculated  x-,  y-,  and  z-  component  of  the  ExB  drift  velocities  near  the  North  Pole. 
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Figure  17:  Trajectories  of  particles  in  the  magnetic  field  acting  under  ExB  drift  with  different  starting  locations  near  the  North  Pole. 
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3.2  Data  Assimilation  Results 


We  performed  a  number  of  data  assimilation  runs  for  different  days  and  months.  Here  we  present 
several  typical  visualization  results  and  later  describe  various  statistical  diagnostics  and  system 
performance  indicators. 

Figure  1 8  and  Figure  19  are  similar  to  the  results  of  pure  model  integration  shown  in  Figure  10  and 
Figure  1 1 .  Note  however  that  these  plots  show  noticeable  local  irregularities  in  response  to  data  inputs. 


Figure  18:  Examples  of  3-D  instantaneous  iso-surfaces  of  electron  density  after  assimilation. 
The  radial  direction  is  exaggerated  to  show  detail. 
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Figure  19:  Examples  of  total  vertical  electron  content  after  assimilation. 


Time  series  of  slant  TECs  from  the  assimilation  and  actual  measurements  are  shown  in  the  next  figure 
for  six  different  reference  stations  used  in  the  assimilation.  These  plots  are  typical  for  all  the  stations  we 
used  and  as  one  can  see  they  demonstrate  agreement  to  within  a  couple  of  TEC  units. 
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Figure  20:  Typical  time  series  of  slant  TEC  from  the  assimilation  for  all  satellites  in  view  (different  colors)  for  six  different  reference  stations. 
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To  test  the  assimilation  system  capabilities  we  withhold  slant  TEC  measurements  for  one  or  more  of 
the  IGS  reference  stations  from  assimilation.  As  data  from  other  stations  is  being  ingested  into  the 
model  we  sample  the  three-dimensional  model  files  of  electron  densities  along  the  line  of  sight  between 
the  control  station  and  all  satellites  in  view. 

We  only  consider  satellites  with  elevation  angles  higher  then  1 5  degrees.  Slant  TECs  computed  this 
way  are  plotted  for  one  station  and  several  satellites  in  view  as  a  function  of  time  in  Figure  2 1 .  The 
assimilation  was  turned  on  at  the  beginning  of  the  integration  and  shortly  afterwards  the  agreement 
between  measurements  and  data  improves.  Depending  on  time  of  day,  satellite  elevation  angle,  and 
number  of  visible  satellites  the  agreement  between  our  nowcasts  and  control  data  proves  to  be  as  good 
as  1-2  TEC  units. 


Figure  21 :  An  example  of  time  evolution  of  slant  TEC  from  the  assimilation  system  (solid  lines)  and  GPS  reference  station 
measurements  (dotted  line)  to  several  GPS  satellites  in  view. 

Note  improving  agreement  after  assimilation  is  turned  on  at  the  beginning  of  the  time  series. 


An  important  diagnostic  of  the  correctness  of  the  assimilation  system  is  the  so  called  7/  (chi-squared) 
test.  The  test,  described,  for  instance,  in  Khattatov  et  al  (2000),  can  be  thought  of  as  a  ratio  between 
Kalman  filter-estimated  background  (model)  error  and  RMS  residual  between  model  and  observations. 
As  such,  if  tunable  parameters  in  the  assimilation  scheme  are  chosen  properly,  this  ratio  should  be  on 
average  close  to  1  after  initial  assimilation  spin-up  time. 

Figure  22,  below,  shows  an  example  of  time  series.  As  one  can  see,  average  y2  is  indeed  close  to  1 . 
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Figure  23:  Average  absolute  root-mean  squared  (RMS)  error  for  assimilation  experiments  with  two  different  sets  of  de-correlation  lengths. 
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To  evaluate  the  system  we  generated  several  measures  of  the  assimilation  system  performance  for  all 
available  reference  stations: 

-  Average  relative  root-mean-squared  differences  (as  a  fraction)  between  simulated  slant  TEC  for 
each  receiver-satellite  pair  and  unbiased  slant  TEC  measurement.  This  diagnostic  is  shown  in  Figure 
24. 

-  Average  absolute  root-mean-squared  differences  in  TEC  units  between  simulated  slant  TEC  for  each 
receiver-satellite  pair  and  unbiased  slant  TEC  measurement.  This  diagnostic  is  shown  in 

Figure  23. 

-  Maximum  or  worst  absolute  differences  between  simulated  slant  TEC  for  each  receiver-satellite 
pair  and  unbiased  slant  TEC  measurement.  This  diagnostic  is  shown  in  Figure  25. 


Average  Relative  RMS  Error,  Lp=0.1*exp(p/10);  Lq=0.01‘exp(p/10);  LLon=0.6 


Figure  24:  Average  relative  root-mean-squared  (RMS)  error  for  assimilation  experiments  with  two  different  sets  of  de-conelation 

lengths. 

In  the  process  of  finding  a  good  set  of  tunable  parameters  for  the  assimilation  scheme  we  performed 
assimilation  experiments  with  several  different  sets  of  de-correlation  lengths.  The  effect  of  different 
parameters  on  system  performance  is  clearly  seen  in  the  abovementioned  plots.  While  we  plan  to 
perform  a  more  systematic  search  of  the  parameter  space  later,  so  far  it  appears  that  the  values 
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corresponding  to  the  second  set  of  plots  are  close  to  the  optimal.  They  result  in  average  RMS  error  of 
between  1  and  2  TEC  units,  relative  RMS  error  of  about  7%,  and  maximum  RMS  residual  of 
approximately  5  TEC  units. 

To  conclude  this  section  we  note  that  the  quality  of  the  produced  electron  densities  depends  critically 
on  the  number  of  reference  stations  used  in  the  assimilation.  We  currently  use  a  list  of  147 IGS 
reference  stations  around  the  world.  Quite  often  their  data  comes  too  late  or  their  RINEX  files  are 
corrupt.  Therefore,  normally  not  all  these  stations  are  used  in  the  assimilation.  A  typical  plot  depicting 
locations  of  a  sub-set  of  these  147  IGS  stations  that  was  used  in  the  assimilation  process  on  a  particular 
day  is  shown  in  Figure  26. 


Worst  (Maximum)  Absolute  Error,  TEC  units;  Lp=0.1*exp(p/10);  Lq=0.01'exp(p/10);  LLon=0.6 
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Figure  25:  Maximum  (worst)  absolute  error  for  assimilation  experiments  with  two  different  sets  of  de-correlation  lengths. 


36 


Figure  26:  Locations  of  IGS  Stations  Used  in  the  Assimilation  Process. 


3.3  Ensemble  Forecasting  and  Targeted  Observations 

We  carried  out  the  ensemble  forecasting  with  various  size  of  ensembles  (5,  1 0  and  20  members)  and 
various  strength  of  the  initial  perturbations  (5%,  10%,  and  20%)  of  electron  densities.  The  control  run 
corresponds  to  the  model  integration  from  March  4, 2004  (0  UT)  to  March  5, 2004  (0  UT).  Every  hour 
we  applied  the  renormalization  of  the  perturbations  for  each  ensemble  member.  Figure  27  and  Figure 
28  present  distributions  of  TEC  fields  after  2  and  6  hour  integrations  and  distributions  of  the  singular 
values  of  the  bred-vector  covariance  matrices  built  with  a  1 0-member  ensemble  size  with  various 
strength  of  initial  perturbations  (20%,  1 0%,  and  5%).  The  positions  of  zones  where  the  singular  values 
are  maximized  do  not  depend  on  the  strengths  of  the  initial  perturbations.  These  positions  typically 
correspond  to  the  sharp  spatial  gradients  of  TEC  and  maximum  ExB-drift  velocities.  “Noisy”  areas  on 
these  plots  are  due  to  small  ensemble  size  dictated  by  computational  requirements. 
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Figure  27:  The  distributions  of  TEC  (a),  leading  singular  vectors  of  C-matrix  for  20  %  (b),  10%  (c),  and  5  %  (d)  for  initial  stochastic 
perturbations  of  electron  density.  Results  correspond  to  the  10  member  ensemble  integrations  after  2  hours. 
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Figure  28:  The  distributions  of  TEC 
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Figure  29  illustrates  the  longitude-altitude  structure  of  the  electron  densities  and  singular  vectors  of  the 
C-matrix  along  the  geographical  equator.  We  can  see  that  the  shape  of  the  singular-vector  distributions 
does  not  depend  on  the  strength  of  the  initial  stochastic  perturbations.  The  maximum  singular  values 
are  located  inside  the  zone  of  the  maximum  vertical  gradients  of  concentrations.  The  transport  of  the 
perturbations  due  to  ExB  advection  plays  a  major  role  in  the  spreading  stochastic  oscillations  in  vertical 
direction.  Figure  30  shows  the  altitude-latitude  structures  of  the  electron  densities  and  corresponding 
singular  vectors  after  6-hour  simulations.  Again  the  spatial  structure  of  the  singular  vectors  depends 
weakly  on  the  strength  of  the  initial  perturbations. 

It  is  worthy  to  note  that  our  results  depend  weakly  on  the  ensemble  size.  For  ensembles  with  5  and  20 
members  (not  shown)  we  have  a  similar  distribution  of  the  leading  singular  vectors  calculated  from  the 
covariance  matrices.  This  leads  us  to  preliminary  conclusion  that  the  ensemble  perturbation  forecasting 
with  the  finite  number  of  members  (10-20)  can  be  employed  to  diagnose  the  regions  for  the  adaptive 
observations  of  electron  density  in  the  ionosphere. 
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Figure  29:  The  longitude-altitude  structure  of  the  electron  density  (a),  leading  singular  vectors  of  C-matrix  for  20  %  (b),  10%  (c),  and 
5  %  (d),  initial  stochastic  perturbations  of  electron  density.  Results  correspond  to  the  10  member  ensemble  integrations  after  2 

hours  at  the  geographical  equator. 

The  nonlinear  plasma  continuity  equations  that  employ  to  forecast  electron  density  distribution  are 
solved  along  the  magnetic  tube  coordinate  system.  From  a  practical  standpoint  we  decided  to  simply 
start  with  analysis  of  the  geographically  regridded  fields.  This  would  be  more  convenient  for 
determining  space  volumes  and  time  windows  for  where  and  when  the  electron  density  sensors  should 

fly- 
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In  the  recent  studies  the  TEC  data  from  about  90  GPS  receivers  provide  a  main  data  stream  for  data 
analysis  and  forecast  of  the  electron  density  fields  in  the  ionosphere.  Daily  about  1 00,  000-200,000  data 
points  (with  time  window  every  5  minutes)  can  be  used  for  global  assimilation  of  the  GPS  TEC 
retrievals.  The  practical  strategy  for  the  adaptive  insertion  of  additional  GPS  receivers  can  be  invented 
in  the  course  of  the  assimilation  of  existing  retrievals.  Depending  on  the  data  analysis  scheme 
(variational  or  sequential)  several  strategies  can  be  proposed  to  evaluate  the  influence  function  of  the 
inserted  'adaptive’  receivers  on  the  initial  forecast  state.  In  this  initial  study  we  decided  to  build  upon  on 
the  ensemble  framework  to  begin  study  of  forecast  errors  using  the  breeding  approach. 

Singular  values  of  the  forecast  error  covariance  matrices  based  on  the  bred  vector  distribution  have 
been  used  to  trace  out  the  most  ‘erroneous’  zones  for  a  given  one  hour  time  window.  Despite  on  the 
rather  complex  and  time  consuming  analysis  of  the  ensemble  simulations  (their  dependence  on  the 
‘starting’  time  of  integration)  the  simple  interpretation  of  results  can  be  proposed  for  given  forecast 
drivers.  These  drivers  are  the  EUV  solar  fluxes,  neutral  wind,  temperature,  and  composition,  and 
empirical  ExB  magnetic  tube  drift  model.  In  particular,  in  the  equatorial  F-region,  a  vertical  ExB  drift 
velocities  in  the  zones  of  strong  vertical  gradients  of  the  plasma  concentrations  create  a  vertically 
elongated  zones  of  highly  uncertain  forecast  predictions.  Near  the  equator,  at  sunset,  the  positions  of 
these  zones  coincide  with  the  locations  of  regions  of  the  Rayleigh-Taylor  instability  ( Kelley  and 
Maruyama,  1 992).  We  can  expect  that  for  the  coupled  ionosphere-thermosphere-electrodynamics 
model  with  the  nonlinear  feedbacks  the  proposed  ensemble  breeding  strategy  can  provide  a  useful 
approach  for  predictability  of  the  equatorial  F-spreads  as  well  as  can  help  to  design  the  special 
‘viewing’  strategies  for  adaptive  observations  during  the  sunset  time  windows. 

In  the  current  uncoupled  formulation,  the  forecasted  plasma  concentrations  depend  on  the  selected 
specification  of  the  neutral  thermodynamics  and  composition  of  the  thermosphere  as  well  as 
distribution  of  the  electrical  and  magnetic  fields.  Specifications  of  these  input  parameters  can 
dramatically  affect  the  ensemble  forecasting  and  change  the  positions  of  the  most  uncertain  predictions. 
The  sensitivity  analysis  of  the  ensemble  forecasting  to  the  stochastic  perturbations  of  the  major  model 
drivers  can  help  to  evaluate  changes  in  the  position  of  the  most  uncertain  zones.  Although  the  coupled 
forecast  ionosphere-thermosphere-electrodynamics  system  will  be  the  most  appropriate  environment 
for  detection  of  the  unstable  (or  most  erroneous)  zones,  the  ‘driver/forcing’  breeding  with  the 
uncoupled  ionosphere  model  can  be  used  as  preliminary  bridge  to  understand  the  sensitivity  of  the 
plasma  forecast  and  its  error  growth  to  stochastic  variations  in  the  external  sources.  We  plan  to 
continue  this  line  of  research  using  the  ensemble  perturbation  framework  in  the  assimilation  of  the  GPS 
TEC  retrievals. 

As  a  final  remark  we  should  highlight  a  self-consistency  issue  between  the  targeting  method  and 
operational  assimilation  schemes.  For  instance,  using  the  4D-Var  scheme  it  is  logical  to  employ  the 
adjoint-based  targeting  techniques  and  singular  value  algorithms  described  in  Palmer  et  al.  [1998].  In 
the  sequential  ensemble-based  assimilation  schemes  the  utilization  of  the  ensemble  targeting  methods 
such  as  breeding  or  replicative  observations  would  be  the  most  natural  strategies  for  adaptive 
observations  ( Lorenz  and  Emmanuel,  1998).  At  this  preliminary  stage  of  the  adaptive  observation 
strategy  for  the  TEC  retrievals  in  the  ionosphere  there  is  still  much  to  investigate  before  to  suggest 
some  practical  recommendations  for  the  operational  data  assimilation  schemes. 

Along  with  the  model  sensitivity  to  the  initial  state,  the  stochastic  and  deterministic  perturbations  of  the 
major  mechanisms  that  control  ionospheric  plasma  structures  should  be  also  analyzed  in  order  to 
understand  the  relative  impact  of  the  uncertainties  in  the  initial  conditions  and  errors  in  external  drivers 
in  the  TEC  forecast  during  the  day. 
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Figure  30:  The  iatitude-altitude  structure  of  the  electron  density  (a),  leading  singular  vectors  of  C-matrix  for  20  %  (b),  10%  (c),  and  5 
%  (d),  initial  stochastic  perturbations  of  electron  density.  Results  correspond  to  the  10  member  ensemble  integrations  after  2  hours. 


4.  CONCLUSIONS 

This  report  presents  final  results  of  a  2-year  AFRL-sponsored  project  whose  objective  was  to  develop 
advanced  modeling  and  data  assimilation  capabilities  for  the  ionosphere  and  upper  atmosphere.  In  the 
course  of  this  project  Fusion  Numerics  developed  a  new  global  three-dimensional  numerical  model  of 
the  ionosphere,  novel  ionospheric  data  assimilation  software  and  methodology,  and  an  infrastructure  for 
obtaining  and  managing  IGS  reference  station  data  for  the  system. 

We  believe  that  the  designed  system  is  unique  in  its  treatment  of  numerical  implementation  of 
ionospheric  physics,  data  assimilation  implementation  and  a  modem  disciplined  approach  to  software 
engineering. 

To  the  best  of  our  knowledge,  this  is  the  only  ionospheric  modeling  and  assimilation  system  existing  in 
the  world  that  actually  solves  all  momentum,  energy  and  mass  conservation  equations  and  implements 
assimilation  in  model  internal  magnetic  coordinates.  Preliminary  results  indicate  absolute  average 
assimilation  error  between  1  and  2  TEC  units  and  relative  error  of  approximately  7%.  Initial  control 
assimilation  experiments  performed  by  withholding  one  or  more  reference  station  data  from  the 
assimilation  demonstrated  accuracy  within  2  TEC  units. 
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Fusion  Numerics  has  received  additional  funding  to  adapt  the  system  for  assimilation  of  data  from  the 
future  Air  Force  satellite  mission  C/NOFS  (Communication  and  Navigation  Outage  Forecasting 
Satellite).  In  the  course  of  the  new  development  effort  additional  features  will  be  developed  and  added 
into  the  system  and  system  performance  is  expected  to  improve. 

An  operational  prototype  of  the  system  has  been  continuously  working  since  August  2003.  The 
prototype  uses  a  first-principles  three-dimensional  time  dependent  numerical  ionospheric  model  as 
forward  propagator  and  a  computationally  efficient  large-scale  Kalman  filter  for  data  assimilation  and 
bias  estimation.  The  system  automatically  acquires  and  archives  relevant  (near)  real-time  or  delayed 
data  from  IGS  network  stations  and  the  NOAA  SEC  center. 

Differential  code  biases  can  be  estimated  on-line  in  near  real  time  and  are  in  approximate  agreement 
with  IGS  Ionospheric  Working  Group  supplied  biases.  One  important  advantage  of  our  method  is  that 
these  biases  (and  satellite  biases  in  the  future)  can  be  made  available  in  real  time.  Should  unexpected 
changes  happen  at  a  particular  reference  station,  (e.g.,  change  of  receiver  or  antenna),  the  system  should 
be  able  to  identify  such  change  and  automatically  adjust  for  it.  Under  some  circumstances,  however,  the 
bias  determination  module  was  known  to  malfunction,  most  likely  due  to  systematic  model  biases.  We 
are  working  on  adding  model  bias  estimation  to  the  data  assimilation  scheme.  For  the  time  being  we 
recommend  using  DCBs  supplied  by  IGS  operational  centers. 

We  believe  that  these  results  indicate  that  the  pursued  approach  is  viable  and  has  practical  merit.  A 
more  extensive  validation  campaign  is  needed  before  such  system  is  deemed  reliable  and  we  are  in  the 
process  of  setting  up  such  a  campaign  with  several  partners. 

One  important  capability  of  the  system  is  that  of  not  only  nowcasting  but  also  forecasting  ionospheric 
conditions.  While  in  principle  a  forecasting  system  is  easy  to  set  up,  as  it  is  not  that  different  from  the 
current  configuration,  rapidly  changing  ionospheric  conditions  impose  high  computational 
requirements.  In  conventional  weather  prediction  forecasts  are  usually  generated  every  several  hours 
since,  say,  humidity  or  temperatures  do  not  change  significantly  in  between.  Ionospheric  electron 
densities  change  rather  abruptly  at  sunrise  and  sunset.  Quite  often  updates  of  1  - 1 0  minutes  might  be 
required.  Considering  that  each  update  requires  a  separate  forecast,  this  can  result  in  significant 
computational  requirements.  The  problem  of  delayed  data  makes  it  necessary  to  perform  forecasts  fast 
enough  to  cover  several  model  hours  in  a  fraction  of  wall-clock  time.  On  the  other  hand,  these 
computations  are  inherently  parallel  and  can  certainly  be  performed  if  enough  processing  nodes  are 
available. 
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